pacman::p_load(rsample)


boot_percentile_function <- function(alpha, 
                                     boots, 
                                     statistics) {
  alpha <- alpha
  statistics <- statistics
  boots <- boots
  bootstrap_percentiles <- purrr::map(
    .x = boots,
    .f = ~ rsample::int_pctl(
      .data = .x,
      alpha = alpha, 
      statistics = statistics
    )
  )
}


# Define Variables  ####

set.seed(853147)

outcome_vars <- c("news_sharing", 
                  "news_distinguish")
experimental_vars <- c("fastpaced_nudge", 
                       "accuracy_nudge")


# Bootstrap of accuracy and accuracy x Authoritarian Models  ####

tic()
## bootstraps
bootstraps <- data_nested %>%
  ungroup() %>% 
  select(data_type, outcome, data) %>% 
  # subset the data for computational feasibility
  mutate(data = map(data,
                    ~ select(.,
                             any_of(outcome_vars),
                             any_of(authoritarian_vars), 
                             any_of(experimental_vars), 
                             any_of(observational_vars),
                             starts_with("num_"), 
                             matches("index"), 
                             "country_id", 
                             "ID"))) %>%
  # nest data by respondent ID for clustered bootstrap resampling
  mutate(data = map(data, 
                    ~ nest(.x, data = -c("ID")))) %>% 
  #sample 5000 times over each subset of the data
  mutate(boots = map2(
    .x = data,
    .y = outcome, 
    .f = ~ rsample::bootstraps(.x,
                               strata = .x[['country_id']],
                               times = boot_times,
                               apparent = T
    )
  )) %>% 
  # estimate model on each of the splits
  mutate(boots = map2(
    .x = boots,
    .y = outcome, 
    ~ mutate(
      .x,
      # lm accuracy model simple
      lm_accuracy_simple_tidy = pmap(
        list(..1 = .x$splits,
             ..2 = .y), 
        .f = ~ fit_lm_accuracy(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data), 
                       model = "simple", 
                       tidy = TRUE
                     ) 
                   ), 
      # lm accuracy model full
      lm_accuracy_full_tidy = pmap(
        list(..1 = .x$splits,
             ..2 = .y), 
        .f = ~ fit_lm_accuracy(
          outcome = ..2, 
          data = analysis(..1)%>% 
            unnest(data),  
          model = "full", 
          tidy = TRUE
          )
        ), 
      # lm accuracy authoritarian interact simple
      lm_accuracy_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "simple", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars))
                     }
                   )
        ),
      # lm accuracy authoritarian interact full
      lm_accuracy_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "full", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars))
                   }
                )
        ), 
      # ME lm accuracy authoritarian interact simple
      lm_accuracy_ME_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "simple", 
                       tidy = TRUE, 
                       authoritarian_measure = "authoritarian_pca", 
                       me_estimate = TRUE)
        ),
      # ME lm accuracy authoritarian interact full
      lm_accuracy_ME_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "full", 
                       tidy = TRUE, 
                       authoritarian_measure = "authoritarian_pca", 
                       me_estimate = TRUE)
        )      
      )
    )) %>% 
  select(data_type, outcome, boots) %>% 
  mutate(boots = map(boots, 
                     ~ select(., c("splits", "id", matches("tidy"))))
  )


# Confidence Intervals of accuracy Models  ####


# calculate the 90, 95, 99 CIs: 

cis_accuracy <- bootstraps %>%
  select(boots, data_type, outcome) %>%   
  # keep only accuracy models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., -matches("authoritarian")))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_simple = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_simple = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_full = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_full = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>% 
  select(data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
  unnest(boot_pctl)


cis_accuracy %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy.csv")




# Confidence Intervals of accuracy x Authoritarian Models  ####

cis_accuracy_authoritarian <- bootstraps %>%
  select(boots, data_type, outcome) %>%   
  # keep only accuracy interaction models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., splits, id, matches("authoritarian"), -matches("ME")))) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_accuracy_authoritarian_simple = map(boots, 
                                                 ~ unnest_wider(., 
                                                                col = lm_accuracy_authoritarian_simple_tidy) %>% 
                                                   select(-splits, -id, -matches("lm_")) %>% 
                                                   rename_with(.cols = everything(), 
                                                               ~ paste0(.x, "_simple")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_accuracy_authoritarian_simple, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_accuracy_authoritarian_simple_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_accuracy_authoritarian_simple) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_accuracy_authoritarian_full = map(boots, 
                                               ~ unnest_wider(., 
                                                              col = lm_accuracy_authoritarian_full_tidy) %>% 
                                                 select(-splits, -id, -matches("lm_")) %>%
                                                 rename_with(.cols = everything(), 
                                                             ~ paste0(.x, "_full")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_accuracy_authoritarian_full, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_accuracy_authoritarian_full_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_accuracy_authoritarian_full) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  select(data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
    unnest(boot_pctl) %>% 
    mutate(authoritarian_specification = names(boot_pctl)) %>% 
    unnest(boot_pctl)


cis_accuracy_authoritarian %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy_authoritarian.csv")



# Confidence Intervals of accuracy x Authoritarian Models (ME estimates) ####

cis_accuracy_ME_authoritarian <- bootstraps %>%
  select(boots, data_type, outcome) %>%   
  # keep only accuracy interaction models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., splits, id, matches("ME")))) %>%   
  # 90 CI 
  mutate(boot_pctl_90_simple = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
  ) %>%   
    # 95 CI 
    mutate(boot_pctl_95_simple = boot_percentile_function(
      alpha = 0.05, 
      boots = boots, 
      statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
    ) %>%   
    # 99 CI 
    mutate(boot_pctl_99_simple = boot_percentile_function(
      alpha = 0.01, 
      boots = boots, 
      statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
    ) %>% 
    # 90 CI 
    mutate(boot_pctl_90_full = boot_percentile_function(
      alpha = 0.10, 
      boots = boots, 
      statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
    ) %>%   
    # 95 CI 
    mutate(boot_pctl_95_full = boot_percentile_function(
      alpha = 0.05, 
      boots = boots, 
      statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
    ) %>%   
    # 99 CI 
    mutate(boot_pctl_99_full = boot_percentile_function(
      alpha = 0.01, 
      boots = boots, 
      statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
    ) %>% 
    select(data_type, outcome, starts_with("boot_pctl")) %>% 
    pivot_longer(cols = starts_with("boot_pctl"), 
                 names_to = c("ci_level", "model_specification"), 
                 names_pattern = "boot_pctl_(.*)_(.*)",
                 values_to = "boot_pctl") %>% 
    unnest(boot_pctl)


cis_accuracy_ME_authoritarian %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy_ME_authoritarian.csv")


toc()







# By Country: Bootstrap of accuracy and accuracy x Authoritarian Models  ####
rm(bootstraps, 
   cis_accuracy, 
   cis_accuracy_authoritarian,
   cis_accuracy_ME_authoritarian)



tic()
## bootstraps
bootstraps <- data_nested_by_country %>%
  ungroup() %>% 
  select(country, data_type, outcome, data) %>% 
  # subset the data for computational feasibility
  mutate(data = map(data,
                    ~ select(.,
                             any_of(outcome_vars),
                             any_of(authoritarian_vars), 
                             any_of(experimental_vars), 
                             any_of(observational_vars),
                             starts_with("num_"),
                             matches("index"), 
                             "country_id", 
                             "ID"))) %>%
  # nest data by respondent ID for clustered bootstrap resampling
  mutate(data = map(data, 
                    ~ nest(.x, data = -c("ID")))) %>% 
  #sample 5000 times over each subset of the data
  mutate(boots = map2(
    .x = data,
    .y = outcome, 
    .f = ~ rsample::bootstraps(.x,
                               times = boot_times,
                               apparent = T
                               )
  )) %>% 
  # estimate model on each of the splits
  mutate(boots = map2(
    .x = boots,
    .y = outcome, 
    ~ mutate(
      .x,
      # lm accuracy model simple
      lm_accuracy_simple_tidy = pmap(
        list(..1 = .x$splits,
             ..2 = .y), 
        .f = ~ fit_lm_accuracy(
          outcome = ..2, 
          data = analysis(..1) %>% 
            unnest(data),  
          model = "simple", 
          tidy = TRUE, 
          by_country = TRUE
        ) 
      ), 
      # lm accuracy model full
      lm_accuracy_full_tidy = pmap(
        list(..1 = .x$splits,
             ..2 = .y), 
        .f = ~ fit_lm_accuracy(
          outcome = ..2, 
          data = analysis(..1) %>% 
            unnest(data),  
          model = "full", 
          tidy = TRUE, 
          by_country = TRUE
        )
      ), 
      # lm accuracy authoritarian interact simple
      lm_accuracy_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "simple", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars), 
                       by_country = TRUE)
                   }
        )
      ),
      # lm accuracy authoritarian interact full
      lm_accuracy_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "full", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars), 
                       by_country = TRUE)
                   }
        )
      ),
      # ME lm accuracy authoritarian interact simple
      lm_accuracy_ME_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ fit_lm_accuracy_authoritarian(
          outcome = ..2, 
          data = analysis(..1) %>% 
            unnest(data),  
          model = "simple", 
          tidy = TRUE, 
          authoritarian_measure = "authoritarian_pca", 
          by_country = TRUE, 
          me_estimate = TRUE)
      ),
      # ME lm accuracy authoritarian interact full
      lm_accuracy_ME_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ fit_lm_accuracy_authoritarian(
          outcome = ..2, 
          data = analysis(..1) %>% 
            unnest(data),  
          model = "full", 
          tidy = TRUE, 
          authoritarian_measure = "authoritarian_pca", 
          by_country = TRUE, 
          me_estimate = TRUE)
      )     
    )
  )) %>% 
  select(country, data_type, outcome, boots) %>% 
  mutate(boots = map(boots, 
                     ~ select(., c("splits", "id", matches("tidy"))))
  )


# By Country: Confidence Intervals of accuracy Models ####

cis_accuracy_by_country <- bootstraps %>%
  select(country, boots, data_type, outcome) %>%   
  # keep only accuracy models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., -matches("authoritarian")))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_simple = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_simple = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_simple_tidy")
  ) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_full = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_full = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_full_tidy")
  ) %>% 
  select(country, data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
  unnest(boot_pctl)


cis_accuracy_by_country %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy_by_country.csv")



# By Country: Confidence Intervals of accuracy x Authoritarian Models ####


cis_accuracy_authoritarian_by_country <- bootstraps %>%
  select(country, boots, data_type, outcome) %>%   
  # keep only accuracy interaction models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., splits, id, matches("authoritarian")))) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_accuracy_authoritarian_simple = map(boots, 
                                                 ~ unnest_wider(., 
                                                                col = lm_accuracy_authoritarian_simple_tidy) %>% 
                                                   select(-splits, -id, -matches("lm_")) %>% 
                                                   rename_with(.cols = everything(), 
                                                               ~ paste0(.x, "_simple")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_accuracy_authoritarian_simple, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_accuracy_authoritarian_simple_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_accuracy_authoritarian_simple) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_accuracy_authoritarian_full = map(boots, 
                                               ~ unnest_wider(., 
                                                              col = lm_accuracy_authoritarian_full_tidy) %>% 
                                                 select(-splits, -id, -matches("lm_")) %>%
                                                 rename_with(.cols = everything(), 
                                                             ~ paste0(.x, "_full")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_accuracy_authoritarian_full, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_accuracy_authoritarian_full_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_accuracy_authoritarian_full) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  select(country, data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
  unnest(boot_pctl) %>% 
  mutate(authoritarian_specification = names(boot_pctl)) %>% 
  unnest(boot_pctl)


cis_accuracy_authoritarian_by_country %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy_authoritarian_by_country.csv")



# By Country: Confidence Intervals of accuracy x Authoritarian Models (ME estimates) ####

cis_accuracy_ME_authoritarian_by_country <- bootstraps %>%
  select(country, boots, data_type, outcome) %>%   
  # keep only accuracy interaction models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., splits, id, matches("ME")))) %>%   
  # 90 CI 
  mutate(boot_pctl_90_simple = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_simple = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_simple = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_simple_tidy")
  ) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full = boot_percentile_function(
    alpha = 0.10, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
  ) %>%   
  # 95 CI 
  mutate(boot_pctl_95_full = boot_percentile_function(
    alpha = 0.05, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
  ) %>%   
  # 99 CI 
  mutate(boot_pctl_99_full = boot_percentile_function(
    alpha = 0.01, 
    boots = boots, 
    statistics =  "lm_accuracy_ME_authoritarian_full_tidy")
  ) %>% 
  select(country, data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
  unnest(boot_pctl)


cis_accuracy_ME_authoritarian_by_country %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_accuracy_ME_authoritarian_by_country.csv")


toc()

rm(bootstraps, 
   cis_accuracy_by_country, 
   cis_accuracy_authoritarian_by_country,
   cis_accuracy_ME_authoritarian_by_country)

